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MRI GRADIENT WAVEFORM DESIGN USING CONVEX 

OPTIMIZATION 

STATEMENT AS TO RIGHTS TO INVENTIONS MADE UNDER FEDERALLY 
SPONSORED RESEARCH OR DEVELOPMENT 

[0001] The U.S. Government has rights in the disclosed invention pursuant to NIH Grant 
Nos. NIH.HL39297, NIH-HL56394, NrH-AR46904 and NIH-CA50948 to Stanford 
University. 

BACKGROUND OF THE INVENTION 

[0002] This invention relates generally to magnetic resonance imaging (MRI), and more 
particularly, the invention relates to the design of magnetic gradients for use in MRI. 

[0003] Magnetic resonance imaging (MRI) requires placing an object to be imaged in a 
static magnetic field, exciting nuclear spins in the object within the magnetic field, and then 
detecting signals emitted by the excited spins as they precess within the magnetic field. 
Through use of magnetic gradient and phase encoding of the excited magnetization, detected 
signals can be spatially localized in three dimensions. 

[0004] Advances in MR system hardware enable the use of new rapid pulse sequences. 
Typical rapid sequences include rapid gradient-echo (FLASH, GRASS) and multi-echo spin- 
echo (TSE, FSE, RARE) sequences. Refocused SSFP (True-FISP, FIESTA, balanced-FFE) 
sequences, which produce high signal and contrast, are becoming common as improved 
gradients allow imaging with minimal artifacts from off-resonance. All of these sequences 
demand efficient gradient waveform design. Efficient acquisition methods include echo- 
planar and spiral imaging. Aside from imaging trajectories, gradient waveform design 
includes phase-encoding, prewinder and rewinder gradients. In rapid sequences with short 
repetition times, the design of these latter gradients is an important consideration in 
improving imaging efficiency, because their duration reduces the image acquisition duty 
cycle. 

[0005] In particular, the design problem is to minimize gradient waveform durations 
subject to both hardware constraints and sequence constraints such as desired gradient area. 
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Numerous previous works have presented different methods to optimize gradients in different 
situations. Many of these methods are limited to the design of trapezoidal pulses, and most 
have been demonstrated for ID gradient design. Simonetti et al. presented a general 
technique to minimize ID-gradient properties such as moments, diffusion-sensitivity 2 or 
RMS-current of ID gradients using numerical optimization. See "An Optimal Design 
Method for Magnetic Resonance Imaging Gradient Waveforms," IEEE Trans Med Imaging 
1993; 12:350-360; and "MRI Gradient Waveform Design by Numerical Optimization," Magn 
Reson Med 1993; 29: 498-504. 



BRIEF SUMMARY OF THE INVENTION 

[0006] The design of rapid imaging sequences requires time-optimal magnetic gradient 
waveforms to provide maximum spatial resolution while keeping both repetition times and 
scan times low. In accordance with the invention, constrained optimization is used in 
designing gradients that satisfy gradient-amplitude and slew-rate limitations and additional 
constraints such as beginning and end amplitude, gradient area and higher-order gradient 
moments. 

[0007] The design is formulated as linear equations for a number, N, of discrete-time 
waveforms with a sampling period t. A minimum time gradient requires the identification of 
the smallest value for N for which a design solution exists. A binary search technique is 
employed to minimize the time in identifying the minimum value of Nr which provides a 
solution. 

[0008] The linear equations are then solved using known linear programming, LI -norm 
formulation, or second-order cone programming (SOCP). 

[0009] The design can also include constraints on the gradient waveform to avoid 
overheating and to limit maximum stimulation thresholds in the waveform algorithm. It is 
known that the gradient amplifier, wires, and coils have several heating constraints. Also, the 
time varying magnetic field of a gradient amplifier can induce electric fields of sufficient 
magnitude to excite nerves. These constraints can be addressed in the gradient waveform 
design. 

[0010] The invention and objects and features thereof will be more readily apparent from 
the following detailed description and appended claims when taken with the drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

[0011] Fig. 1 illustrates a voltage model and a constant-slew-rate-limit model. 

[0012] Fig. 2 illustrates discrete approximation of a gradient waveform. 

[0013] Fig. 3 illustrates piece- wise-linear approximation to quadratic constraints. 

[0014] Fig. 4 illustrates piece- wise-linear approximation to a spherical surface. 

[0015] Figs. 5a-5c illustrate comparison of speed of different formulations for ID, 2D and 

3D, respectively. 

[0016] Figs. 6a-6d illustrate first-moment-nulled gradient design. 

[0017] Figs. 7a-7b illustrate echo-planar transition waveform design for gradient 

waveforms and slew-rates. 

[0018] Figs. 8a-8d illustrate spiral rewinder gradient design. 
[0019] Figs. 9a-9d illustrate moment-nulled spiral rewinder gradient design. 
[0020] Figs. 10a- 10c illustrate "stack-of-spiral" rewinder gradient design. 
[0021] Fig. 1 1 illustrates gradient power for different oblique gradient angles. 
[0022] Figs. 12a- 12b illustrate on-axis echo-planar transition waveform design. 
[0023] Figs. 1 3a- 1 3b illustrate off-axis echo-planar transition waveform design. 
[0024] Fig. 14 is a graph illustrating voltage and current constraints and magnetic 
stimulation constraints for safe operation of a gradient amplifier. 



DETAILED DESCRIPTION OF ILLUSTRATIVE EMBODIMENTS 

[0025] The design method can be used in accordance with the invention to design 
minimum-time gradient waveforms that satisfy gradient-amplitude and slew-rate limitations. 
Additional constraints such as gradient start and end amplitudes, gradient area or higher-order 
gradient moments can be included. Importantly, application of the method can be extended 
to multi-dimensional gradient design which is necessary for procedures that use oblique scan- 
plane orientations. 

[0026] Consider first the general constraints in gradient waveform design including 
gradient amplifier voltage and current limits as well as pulse-sequence requirements. Next, 
methods of formulating the problem are described so that it can be solved using three 
different types of existing linear and quadratic programming techniques. All of these 
techniques result in comparable solutions. 
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[00271 A gradient amplifier supplies current to a gradient coil, producing a continuously- 
varying gradient, G {t). The amplifier has current and voltage limits (W and Vmax) that 
result in the following limits on G (0 as described in King et al., "Optimized Waveforms for 
Spiral Scanning," Magn Reson Med 1995; 34:156-160. 



and 



L — G(t) + RG(t) 
dt 



< nV 

— max 



[0028] Here L,Rmd^ are the gradient coil inductance, resistance and efficiency (i.e., in 
mT/m/A). hi MRI systems, the gradient can be considered a vector, G(0 that is composed 
of three components, G, (t), Gy (t), and G, (t). The primary purpose of the invention is to 
specify gradient waveforms that can be freely rotated in three dimensions. This constrains 
the maximum gradient or change of gradient in any direction, and the constraints of Eqs. 1 
and 2 become (with Gmax = n^max) 



\G(,T)1^<iiG„ 



and 



(3) 



iL^-GiO + RGiti <TjV^ (4)- 
II dt II 2 

[0029] Frequently, a simpler "constant-slew-rate-limit" model is used, where r\RImax is 
subtracted fi-om the right side of Eqs. 2 and 4 and the RG (t) or R Git) term is omitted. For a 
single gradient axis, the constraints of the voltage model and the constant-slew-rate-limit 
model are shown in Fig. 1 which illustrates limits on gradient amplitude and slew-rate for a 
single-axis gradient due to maximum current (Imax) and maximum voltage (Vmax). Limits are 
shown using the common constant-slew-rate-limit model (white area) and the more accurate 
voltage-limit model (white plus light gray areas). As the ratio of coil resistance (R) to coil 
inductance (L) increases, the constant-slew-rate-model becomes increasingly constraining, 
and the voltage-limit model should be used, rj is the ratio of gradient strength to gradient 
current. As the ratio of coil resistance to coil inductance increases, the constant-slew-rate- 
limit model becomes more constraining. 

[0030] Imaging requirements of the pulse sequence constrain the k-space change ( M ) that 
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the gradients will produce in terms of the gyromagnetic ratio, y as 

Zo^i^y^ (5)- 

Additionally, there are cases where higher-order moments of the gradients can be specified, 



as 



ly^i^y' =-^^^0 (6) 

where Am^ is a specified design constraint on the change in the qth moment of the 
gradient. 

[0031] Boundary constraints on the gradient waveforms include having specified initial and 

final values for the gradient, i.e., G(o) = and G{t) = gj. . This is required to design 

gradient waveforms at the start and end of the sequence, and to link different waveforms 
together. 

[0032] Practical gradients are generated as discrete-time waveforms with a sampling period 
t. A gradient waveform can be represented as a discrete sequence G[n] = G(«r) with n = 

Figure 2 shows a single component of the discrete-time vector G[n] . The constraints 
in Eqs. 1-6 can easily be converted to discrete-time equivalents, where a derivative is 
approximated by the slope between two consecutive samples, and an integral is approximated 
by a discrete sum multiplied by x. Specifically, it is usefiil to express the voltage constraint 
(Eq. 4) as 

aG[n]+ j8G[n-hl]|| <^ (7) 
lb 

where a=i?- — ,)9 = -^, and 0 = rjV_^^ for the voltage model. In the constant-slew- 
r r 

rate-limit model, a = -^,y5 = -,andi/^ = rjV^ - RI^^ . 

r T 

[0033] Typically the gradients, as well as the desired moments, consist of two or three 
separate gradient axes, hi basic cases, these axes can be treated separately and all constraints 
expressed as described in the previous section. However, there are numerous cases where 
gradients must be designed in a 2D or 3D space that must be fi-ee to rotate. 

[0034] First we now define the solution space consisting of three separate discrete-time 
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waveforms along each axis, [n], Gy [n] and G, [n], which must collectively satisfy the 
quadratic constraints of Eqs. 3 and 4. The boundary constraints and moment constraints 
(Eqs. 5 and 6) are applied separately for each of the dimensions. Overall, the problem is no> 
quadratically constrained. 

[0035] Before formulating solutions to the gradient design problem, it is useful to 
summarize the constraints for discrete time waveform design. Table 1 lists the four types oi 
constraint for gradient design problems for D-dimensional, N-point waveform with Q- 
moments constramed, and the number of instances of each in gradient design. The quadrati 
inequality constraints are the prime source of complexity in gradient design. 

[0036] Table 1. 



Constraint 


Type 


Number 


Current (Gradient Amplitude) 


Quadratic Inequality 


N 


Voltage (Gradient Slew-Rate) 


Quadratic Inequality 


N 


Gradient Moments 


Linear Equality 


DxQ 


Boundary Values 


Linear Equality 


2D 



[00371 A certain class of constrained optimization problems, convex-optimization 
problems, have been studied heretofore. For many types of convex optimization problems, 
very efficient and robust solution methods have been developed. These methods are 
guaranteed to find the globally-optimal solution if any solution exists. Following are 
described methods of expressing the gradient design constraints as standard convex 
optimization problems so that they may be solved efficiently and reliably. 
[00381 Both absolute-value and equality constraints can equivalently be expressed as two 
inequality constraints. A single absolute-value constraint a\x\ < 6 can be replaced by two 
linear constraints, ax < b and - ax < b . An equality constraint of the form ox = can be 
replaced by the two inequalities ax < b + s and - ax < -b + e . This is desirable in practical 
optimization, as it allows the tolerance e to be different for different constraints. 
[00391 Using this formulation, the constraints of Eqs. 1 , 2, 5 and 6 as well as the boundary 
constraints can all be expressed as linear inequality constraints on the discrete-time gradient 
waveform. 
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[0040] Each quadratic constraint may be approximated by multiple linear constraints. This 
is done in 2D, for example, by approximating the circle Gx [n^ + Gy [ri^ = Gm<J as a series of 
P lines defined by 



cos((9^)G,[/i] + sin(e^)G^[«] < G^, cos 

Vp) V = (8) 

where 6^ = In ^ . Figure 3 shows this piece-wise-linear approximation to a 2D 

quadratic constraint with P=8. The white area represents allowable values of the 2D 
amplitude, voltage, or slew-rate. For all of the 2D examples herein, P=16 is used. A piece- 
wise-linear approximation to a sphere is also possible for 3D quadratically constrained 
problems. A simple way to do this is to consider the surface of the sphere | r | = c in the first 
octant where the components of r are all positive (Fig. 4). By placing points on the surface 
of the sphere, triangles with roughly equal area are constructed, and the constraints can be 
expressed in the form 

nr<n'V (9) 

where w is a vector normal to the triangle, i.e., n = (5 - 3) - x(c - ^) and V is the 
position vector of any of the triangle comers, {A, B ,or C). 

[0041] For the 3D examples herein, the first quadrant is divided into 9 planes using the 
points shown in Fig. 4. Here piece-wise-Unear approximation to spherical surface represents 

3D quadratic constraints. The points shown are permutations of (l,0,0) , » 

( ' ^ ' ^) ®^tire sphere can then be represented by the 72 planes generated with 

all permutations of positive and negative signs in the individual n vectors. 
[0042] The problem with only linear inequality constraints can be solved using linear- 
programming (LP). The "standard LP" form of a problem is to find the vector x that 
minimizes a cost function f subject to the constraint Ax < b. The matrix A and vector b 
are formed by combining all of the linear constraint equations for amplifier and pulse 
sequence constraints. This formulation of the problem is referred to as "simple-LP" 
formulation. 
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[0043] The cost function in simple-LP formulation is linear in the gradient values. 
Minimization of such a cost function is not particularly useful. As such the actual 
optimization consists of a series of feasibility tests as will be described later. For the simple- 
LP formulation, f is chosen to be a vector of ones. 

[0044] The LP formulation for an TV-point waveform uses ND variables with roughly 
2 NP constraints, where D is the number of dimensions and P is the number of piece-wise- 
linear constraints on gradient amplitude and voltage amplitude for which coefficients are all 
positive (i.e., the number of constraints in the first quadrant or octant). The Appendix shows 
the complete simple-LP formulation for a 2D freely-rotatable gradient design problem. 

[0045] An alternative to simple-LP formulation is the LI -norm formulation, which alters 
the number of constraints and can produce different solutions to the problem, see Xu et al., 
"Homogeneous Magnet Design Using Linear Programming," IEEE Trans Med Imaging 
2000; 36:476-483. 

[0046] The LI -norm formulation uses standard linear programming to minimize the LI - 
norm (absolute value) of the gradient and slew voltage by adding "slack variables" to the 
optimization problem. In addition to solving for Gx [«], we solve a set of variables Hx [w], 
that converge to \Gx and a set of variables Sx [n] that converge to \aG^ [«]+ fiG^ + 1] • 

[0047] The slack variables are forced to converge by first adding the following linear 
gradient constraints to the problem: 



[0048] Additionally we add the following slew-rate constraints: 

- S,[n\ + aG^[n\ + /3G^{n + 1] < 0 
-5,[«]-aGJn]->ffGJ« + l]<0 

[0049] These constraints will force the Hx [«] and Sx [«] variables to approach the 
appropriate absolute value when combined with minimization of the following cost function: 



[0050] For two and three dimensional problems, an appropriate set of variables Hy [n], Sy 
[«], H2 [n] and Sz [n] would be added. Constraints for these variables are expressed as in Eqs. 



-//Jn] + G,M<0 
-H^[n]-G,[n]<0 



(10). 



(12). 
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10 and 11, and the variables would be added to the cost function in Eq. 12. 

[0051] The simple-LP cost function of all ones simply tries to make individual gradient 
axes closer to their negative amplitude limits. However, the LI -norm cost function is 
preferable, as it tends to minimize the magnitude of both gradient current and voltage, which 
will improve heating characteristics of the waveform. 

[0052] When LI -norm formulation is used, the piece-wise-linear constraints on gradient 
amplitude may be reduced by replacing the constraint on Gx [n], Gy [«], and Gz [n] with 
equivalent constraints on Hx [w], Hy [w], and Hz [n]. Similarly, the voltage constraints may be 
replaced with constraints on Sx [«], Sy [n], and Sz [n]. Since all of the slack variables are 
constrained to be non-negative, only the gradient amplitude and voltage-limit constraints with 
positive coefficients are needed. This reduces the number of gradient amplitude constraints 
by a factor of 2^ where D is the number of gradient axes, and is the primary advantage of the 
LI -norm formulation compared to simple-LP formulation. 

[0053] The LI -norm formulation uses 3ND variables, and approximately 2N(P + 2D) 
constraints. For 2D or 3D problems, the number of constraints is often much smaller than 
that of the simple-LP formulation. Additionally, as described by Xu et al., LI -norm 
programming tends to find solutions where some or many of the variables are zero. This may 
be advantageous in gradient design, as solutions with lower gradient or slew duty cycle result 
in lower gradient heating. 

[0054] Second-order cone programming (SOCP) (see Lobo et al., "Applications of second- 
order cone programming," Linear Algebra Appl. 1998; 284: 193-228) is a method that finds 
the solution x that minimizes a linear cost function f subject to the "second-order cone" 
constraint. 

Il^ + Z^ll^ <Cx'hd (13). 

[0055] First this constraint is a superset of the linear constraints in linear programming. 
Second, the cost function is linear, as with the simple-LP formulation described above, and is 
chosen in the same manner, as all ones. 

[0056] The SOCP formulation of the gradient design problem begins with Eqs. 5-6. 
Although these constraints could be reduced by a factor of 2^ compared with the simple LP or 
LI -norm formulations, this has little effect on convergence time. The main advantage of the 
SOCP formulation is that the quadratic limits of Eqs. 3 and 4 result in a total of IN -l 
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constraints, regardless of the number of gradient dimensions. Additionally, the SOCP 
formulation does not approximate the quadratic gradient limits, which sometimes results in 
slightly shorter duration solutions. 

[0057] Because the time to find a solution increases with the number of constraints, it is 
useful to try to further reduce the number of constraints. Assuming that the slew-rate 
constraints are met, there will be gradient amplitude constraints near the endpoints that are 
redundant. Specifically, if the (minimum) distance from the start or end gradient to a linear 
or quadratic gradient constraint is AG, then that constraint is redundant for at least a time 



[0058] The minimum-time gradient design problem is to find Nmw, the minimum value ofN 
for which a solution exists. Linear programming functions can be called sequentially with 



for N, then Nmin < N, Otherwise, Nmin > N . 

[0059] This suggests using a binary-search technique. This technique divides the unknown 
interval containing N^jin on each call. Begin with estimates of the upper and lower bounds Nu 
and Ni on A^, First these are tested to ensure that they are in fact upper and lower bounds. If 
not, both bounds are raised or lowered accordingly and retested. Once Nu and A^/have been 
verified, A'^ is set as close as possible to Ni + v(7V„ - Ni), Depending on whether or not a 
solution exists, either Nu or Ni is set to A^, and the process is repeated. Use v=0.8 because an 
existing solution is found considerably more quickly than the lack of a solution using SOCP. 

[0060] The above formulation has been implemented using Matlab 6.0 with fiinctions from 
the Mathworks Optimization toolbox (Mathworks, Natick, MA). Specifically, use the 
linprogO function, which replaces the old IpQ function in Matlab, as well as the socpQ 
function described above. Matlab functions are available for general use. 

[0061] First, compare the time required to solve the minimum time gradient-area-nulling 
problem using the simple-LP, LI -norm and SOCP methods. All tests used Matlab 6.0 on L8 
GHz Athlon PCs running Linux. The binary search formulation was identical for all three 
methods, and the comparison was performed for ID, 2D and 3D problems. 

[0062] Generated were 1000 random gradient waveforms of random duration and the given 
dimension with endpoints within the magnitude constraint for the given solution method. A 



max 



different values of A^. We first assume that either G(o) 



0 or G(wr)=0. If a solution exists 
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sample period of r =20 (is was used, which assumes a gradient ampHfier bandwidth of less 
than 25 kHz, The gradient and zeroth moment of the gradient were then rewound to zero 
using the bisection method described above. Each step of the binary search method was 
solved using simple-LP, LI -norm or SOCP formulation. The times required to solve the 
rewinder problem for each waveform and each formulation method were recorded. Li all 
cases a maximxmi gradient amplitude of 40 mT/m, and a constant-slew-rate model with 150 
T/m/s maximum slew-rate were used. For 2D and 3D, the rewinders are limited by the 
quadratic constraints described above. 

[0063] Figures 5a-5c compare convergence times for the three methods for ID, 2D and 3D 
design problems. Shown are convergence time for ID (a), 2D (b), and 3D (c) time-optimal 
gradient design using simple-LP (dotted line), LI -norm (dashed line) and SOCP (solid line). 
Although simple-LP formulation is the most efficient for ID design problems, the LI -norm is 
significantly better for 2D and 3D problems. SOCP outperforms LI -norm in 2D and 3D 
problems, except when the solution tums out to be very long. Not surprisingly, convergence 
time increases with the solution length. SOCP formulation is fastest regardless of dimension 
for most practical waveforms (i.e., duration less than 1 ms), and is significantly faster than 
the other methods as the number of gradient dimensions increases to 2 or 3. 
[0064] Presented now are several examples showing applications of linear programming 
for time optimal gradient design. We begin with ID examples that show that linear 
programming produces the time-optimal waveforms. Next we show 2D examples, 
particularly for use in spiral imaging. Finally we show 3D examples for double-oblique or 
oblique spiral scans. 

[0065] For a constant amplitude limit, constant-slew-rate limit gradient, the minimum time 
waveforms to produce a given zeroth, first and/or second moment are generally triangular or 
trapezoidal waveforms. However, calculation of the segments sometimes must be solved 
numerically. 

[0066] A first example is the generation of a simple ID moment-nulled phase-encoding 
gradient. The boundary constraints on the pulse waveform are Gx [0] = G> [AT] = 0 and the 

desired k-space area and first-moments are A A: = k^^^ and mi = 0. This example sets 

k^^ = 400 m"\ Gmax = 40 mT/m and maximum slew-rate — = 150 T/m/s. 

L 

[0067] The resulting gradient waveform, k-space trajectory, slew-rate and first moment are 
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shown in Fig. 6 which shows (a) first-moment-nulled gradient waveform, (b) corresponding 
k-space trajectory, (c) gradient slew-rate, and (d) gradient first-moment, all plotted as 
fimctions of time for Example A. The gradient waveform has a duration of 1 .048ms. The 
gradient is a bipolar combination of a trapezoidal lobe with a triangular lobe, as can 
reasonably be expected given amplitude- and slew-limited constraints in one dimension. 
Note that either the gradient amplitude limit or the slew-rate limit is always met, which 
should be the case for a minimum-time gradient in one dimension. 

[0068] A second example is to design a minimum-time waveform that transitions between 
successive lines of an echo-planar imaging trajectory. For consistency, the same gradient and 
slew-limits, 20 mT/m and 200 T/m/s respectively, are used with the assumption that the 
gradients can be fi^eely rotated to any oblique in-plane position. The EPI trajectory uses a 20 
mT/m readout gradient, with 0.8 cm*^ phase-encode area between lines. In our formulation, 
we specify [1] = -G^ [TV] = 20 mT/m, Gy [1] = [TV] = 0 mT/m, dikj, = 0 cm"' and A;t^=0.8 
cm \ T = 1 |xs. 

[0069] The result, shown in Fig. 7, has a total duration of 245 |is. In Fig. 7, gradient 
waveforms (a) and slew-rates (b) for an echo-planar readout/phase-encode transition are 
shown. Gradients and slew-rates are shown for the x-axis (solid line) y-axis (dashed line) and 
magnitude (dotted line). This method performs as well as the time-optimal method described 
in ref (9). 

[0070] The solutions with both methods look very similar, and the minor difference 
between them is most likely simply the way that the end-point constraint is formulated. 
Interestingly, this gradient is not zero at either endpoint. However, because of the symmetry 
of the problem, the bisection method of finding the minimum-time gradient still works. 

[0071] Minimum gradient waveform design on one dimension can generally be solved 
exactly and analytically. However, two-dimensional design is not always simple. Consider 
spiral imaging, where interleaved spiral waveforms are played by rotating the waveform in 
the x-y gradient plane. For the trajectory to be played at arbitrary rotation angles, the 
gradient and slew-rate vector magnitudes must be within slew-rate constraints. 

[0072] Spiral waveform design methods have been addressed previously. The spiral 
waveform can be rewound by separately rewinding each axis with both amplitude (current) 
and voltage limits set to 1/ V2 times the individual axis limits. This is adequate for 
applications where the spiral readout duration is long, but is not time-optimal. In cases such 
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as spiral steady-state free precession imaging, it is desirable to rewind the first moment of the 
gradient as well. In this example, a time-optimal rewinder for spiral gradient waveforms as 
well as a rewinder for a first-moment-nulled spiral waveform are used 

[0073] For this example, a spiral waveform can be used similar to that used for SSFP 
imaging. The design assumes 60 spiral interleaves to achieve a resolution of 1 .25 mm and a 
field-of-view of 20 cm with a readout duration of 2.7 ms. 

[0074] The optimization goal is to rewind the gradient waveforms and zeroth moment 
vector simultaneously to zero as quickly as possible. The rewinder gradient vector initially 
equals the final gradient of the readout, and ends at zero in both dimensions. Combined with 
the constraint of rewinding the zeroth moment vector, this gives six equality constraints. 
Figure 8 shows the spiral gradient waveform and rewinder characteristics. Figs. 8a-8d 
illustrate gradient waveforms (a), k-space trajectory (b), slew-rates (c) and first moment (d) 
for a simple rewinder waveform (gray region) for a 2D-spiral imaging gradient. Each plot 
shows X (solid line) and y (dashed line) components, as well as the magnitude (dotted line). 
The k-space-only rewinder lasts for 0.66 ms, and has a significant residual first moment. All 
quantities reach zero at the end of the waveform. As with the ID example, either the gradient 
magnitude or slew-rate magnitude limit is reached for the entire duration of the rewinder. 
[0075] The spiral rewind problem can be repeated for the case where the first moments are 
also rewound at the end of the waveform. This adds a constraint in each axis for a total of 
eight equality constraints. The LI -norm solution is shown in Figs. 9a-9d. Here are shown 
gradient waveforms (a), k-space trajectory (b), slew-rates (c) and first moment (d) for an 
nuUed rewinder waveform (gray region) for a 2D-spiral imaging gradient. Each plot shows x 
(solid line) and y (dashed Hne) components, as well as the magnitude (dotted line). The wi- 
nuUed rewinder lasts for 1.23 ms, and rewinds the gradient, k-space location and the first 
moment. For nulling of the first moments, the minimum-time waveform is significantly 
increased, fi-om 0.78 ms to 1 .38 ms. If the two axes were rewound separately (as is 
commonly performed), each using 70% of the maximum power, the duration would be closer 
to 1 .8 ms. hi rapid imaging sequences such as SSFP, 0.4 ms is a significant time reduction, 
i.e., 10% of the repetition time or 25% of the data acquisition time. 
[0076] There are cases where the gradient vector must be rotated freely in three 
dimensions, such as double-oblique scans or oblique-plane spiral scans. The 3D "stack of 
spirals" example in Figs, 10a- 10c pertains to the latter. Shown are gradient waveforms (a), k- 
space trajectory (b) and slew-rates (c) for a 3D rewinder waveform (gray region) for a "stack- 
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of-spirals" imaging gradient. Each plot shows x (solid line), y (dashed line) and z (dash-dot 
line) components, as well as the magnitude (dotted line). The rewinder gradient lasts for 0.74 
ms. For gradient-spoiled acquisition, it is desirable to play a gradient spoiler on the z-axis 
following the spiral readout. To minimize the overall time, this spoiler should be played out 
simultaneously with the spiral rewinder. 

[0077] The problem is formulated the same way as that of the simple spif al-rewinder in 
Example B with the additional constraints that in the z-axis, we specify Gx[l] = Gx[N] = 0 
and a desired k-space area. 

[0078] The LP formulation is useful for optimizing oblique gradient design. Gradient 
design is simplest when the three gradient axes may be designed independently of each other, 
that is when each "logical" gradient axis corresponds exactly to one "physical" gradient axis. 
This non-oblique case corresponds to the 3D LP formulation where the gradient constraints 
are simply \Gx [n]\ < G„,ax, |G> [n]\ 

— Gmax9 

[0079] Gradient limits can be reset for any rotation angle by scaling the logical gradient 
limit by (cos (p + sin where (p is the angle between major logical and physical axes. As 
shown in Fig. 11, this ensures that the gradient vector always lies within the large gray 
square. However, this is obviously inefficient in that a significant amount of gradient power 
is not used. Fig. 1 1 illustrates gradient amplitude usage for different oblique-scan strategies. 
The dark gray square represents all usable gradient amplitude, and the circle represents the 
quadratically-constrained amplitude if the gradient is designed to be fi-eely-rotatable. If 
gradients are designed completely independently, then the amplitude limits on each axis must 
be reduced so that the worst-case gradient still lies within the dark gray square. This is 
shown for a rotation angle of ^ = 15"* (dashed lines). With constraints tailored to the specific 
oblique angle, all gradient power is available for any rotation angle. 

[0080] The LP formulation is very adaptable, as the physical constraints can simply be 
rotated to the logical coordinate system. This means that for oblique gradient design all 
available gradient power is always used. 

[0081] Consider the echo-planar transition of example B where the logical and physical 
gradients are rotated by ^zJ= 0° and (f>= 30''. Here the design is gradient for each specific 
oblique position. The gradient amplitude constraints are 

COS^,[«]-sin^^[«] < Gmax 
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- cos^j,[n] + sin^^[«] 



< G, 



max 



sin (l>G^ [n] + cos [«] 



< a 



max 



- sin 4<}^ [n] - cos 4<jy [«] 



< a 



max 



(14) 



where ^ is the angle of rotation between logical and physical coordinate systems. The 
slew-rate constraints are 



[0082] Figures 12 and 13 show the resulting waveforms. Both waveforms are faster than 
that shown in Fig. 7, as more gradient power is available, both for amplitude and slew-rate. 
The off-axis case results in a longer waveform, because slewing in two logical axes must be 
performed by one physical axis. 

[0083] The gradient waveform design can also include constraints to avoid overheating and 
to limit maximum thresholds. The gradient amplifier, wires, and coils have several heating 
constraints. These can be expressed as functions of the gradient waveforms and the 
derivatives of the waveforms. Often these constraints are binding in the operation of the 
scaimer. Hence it would be important to be able to constrain the heating in the waveform 
design algorithm. Fortunately it is possible to include constraints on the gradient waveform 
and its derivative in the methods disclosed herein. 

[0084] The heat dissipated in the coil itself is 



where Rcon is a constant, and / is proportional to the gradient waveform, G(t). 

Using the SOCP formulation, the waveforms G(t) for each axis are approximated by a 
discrete-time waveform G/, where /=1 . . .N. Using this, the quadratic constraint can be 
expressed directly and added to the other constraints in the problem, as 



cos <l>{aG^ [«] + y(?G^ [« + 1]) - sin (f>[aGy [n\ + fiG^ [az + 1]) < 
- cos<^(aGJn] + >8GJ« + 1]) + sin ^(aG^[w] + pGy[n + 1]) < 

sin (piaG^ [n] + fiG, [n + 1]) + cos (l>[aGy [«] + fiGy [« + !])< 
-sin<2>(aGJ«] + pG^[n + 1])- cos^(aG^[«] + pGy[n + 1]) < 



(15). 




(16) 



0 
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tZ^Gf <^fnEZ (17) 



1=1 



where Emax is the maximum heat that can be dissipated in the coil for a given time duration, r| 
is the coil efficiency (G/cm/A), and Tis the gradient sampling step. Alternatively, instead of 
constraining the heat, the heat can be minimized by first adding a slack variable, E, in the 
SOCP formulation as follows: 



\ZtG]<E (18). 



1=1 



Then E can be minimized in the SOCP cost function directly. 

[0085] Amplifier heating models vary firom system to system, but generally are of the 
following form: 

T T 

Ea.p = \l"R.a.dt+ \Wsa,dt+V^Jl\PWM{I) (19). 
0 0 

The first two terms of this equation can be constrained or minimized by adding slack 
variables in the SOCP formulation. The first term uses one slack variable, Ej, as 



2 



2LtR,„,^<E, (20). 
The second term first defines variables, Pi, as 



'^V-<^- (21). 



Then, neglecting the third term, Eamp can be constrained or minimized by constraining or 
minimizing the sum: 

N 

J = '£tP,+E, (22). 

1=1 

[0086] The time varying magnetic field of a gradient amplifier has been known to induce 
electric fields of sufficient magnitude to excite nerves. This is known as magnetostimulation 
or peripheral nerve stimulation (PNS). Stimulation is a fimction of gradient amplitude and 
the derivative of the gradient amplitude. The FDA does not advise MRI manufacturers to 
operate the gradient amplifiers above stimulation thresholds. Hence, it would be a 
contribution to the state of the art of gradient waveform design to be able to include 
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maximum stimulation thresholds in the wavefomi algorithm. Average population studies 
have determined stimulation thresholds that are a linear combination of G(t) and dG/dt(t). 

[0087] Mathematically, this stimulation limit can be written as 



aG{t)^P- 



dt 



<Astin. (23) 



where G(t) is the gradient and dO/dt is the first derivative of the gradient. The stimulation 
threshold, Asum^ is determined experimentally, and may vary from patient to patient. Given 
Astim-, the stimulation limits can be imposed on the individual or combined gradient axes in the 
exact same way as the amplifier voltage limits. Fig. 14 illustrates the voltage and current 
constraints and magnetic stimulation constraints for safe operation of a gradient amplifier. 

[0088] The design of many minimum-time gradient wavefomis can be expressed as a 
constrained optimization problem. With some manipulation, the problem can be posed as a 
standard linear programming (LP) problem, for which many efficient solving methods have 
been developed. The formulation described solves minimum-time gradient problems where 
the initial and final gradient values and desired moments are specified, and where the 
gradient either starts or ends at zero. 

[0089] Proof of time-optimality is an interesting question. For ID cases, this can be done, 
and results indicate this algorithm produces the minimum-time solution. For more complex 
cases where there is not an analytic solution, it is not possible to prove time-optimality of the 
solution. Based on the principles of convex optimization techniques, it is assumed that the 
solutions are time-optimal. Additionally, the simple LP formulation and LI -norm 
formulation always give the same duration gradient. 

[0090] The minimum-time gradient problem, as expressed above, is a series of feasibility 
problems using different length solution vectors. Shown are three different problem 
formulations. First, the problem can be solved simply as a standard LP problem, using 
\n\ Gy [«], and G^ \n\ variables representing the gradient on each axis. Second, the LI -norm 
formulation adds additional slack variables, Hj, [n], Hy [w], and [«], representing the 
gradient magnitudes, and Sx [«], Sy [«], and Sz [«] representing the voltage magnitude to 
fiirther reduce the number of constraints. The LI -norm formulation also allows a more 
physically-meaningfiil cost fimction to be expressed. Third, the second-order cone 
programming (SOCP) formulation uses Gx [«], Gy [«], and G^ [«] variables and quadratic 
constraints on both gradient amplitude and gradient voltage. 
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[0091] The simple LP formulation is advantageous compared to Ll-norm only for ID 
problems. For 2D or 3D problems for which the solutions consist of 75 or fewer time 
samples (1.5 ms with Ar=20 ^is), the SOCP formulation is the fastest formulation. Given that 
the SOCP formulation also uses the available gradient limits slightly more efficiently than 
either simple LP or Ll-norm formulation, it is clearly the best option. However, if SOCP 
software is not available, the Ll-norm formulation also performs very well for 2D or 3D 
problems. 

[0092] Minimum time, multi-dimensional gradient design is an increasingly important part 
of designing MR imaging sequences. The invention provides a method that can deliver the 
time-optimal gradients given a variety of constraints such as amplitude, slew-rate, end points 
or gradient moments. The invention uses standard linear programming tools, which are 
robust and efficient. 

[0093] While the invention has been described with reference to specific embodiments, the 
description is illustrative of the invention and is not to be construed as limiting the invention. 
Various modifications and applications may occur to those skilled in the art without departing 
fi-om the true spirit and scope of the invention as defined by the appended claims. 

[0094] While the invention has been described with reference to specific embodiments, the 
description is illustrative of the invention and is not to be construed as limiting the invention. 
Various modifications and applications may occur to those skilled in the art without departing 
fi-om the tme spirit and scope of the invention as defined by the appended claims. 
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